{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mesh generation for tutorial 06, case 7b\n",
    "\n",
    "This file generates the mesh which is used in the following examples:\n",
    "* 7b_stokes_neumann_control\n",
    "\n",
    "The test case is from section 5.5 of\n",
    "F. Negri. Reduced basis method for parametrized optimal control problems governed by PDEs. Master thesis, Politecnico di Milano, 2010-2011."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import gmsh\n",
    "from mpi4py import MPI\n",
    "from dolfinx.io import XDMFFile\n",
    "from multiphenicsx.mesh.utils import gmsh_to_fenicsx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert MPI.COMM_WORLD.size == 1, \"This mesh generation notebook is supposed to be run in serial\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometrical parameters and related quantities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1 = 1.0\n",
    "mu2 = np.pi / 5.0\n",
    "mu3 = np.pi / 6.0\n",
    "mu4 = 1.0\n",
    "mu5 = 1.7\n",
    "mu6 = 2.2\n",
    "lcar = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = 1.0\n",
    "X = -Y\n",
    "L = 3.0\n",
    "B = Y - mu1\n",
    "H_1 = B + np.tan(mu2) * mu5\n",
    "H_2 = B - np.tan(mu3) * mu6\n",
    "L_1 = mu1 * np.cos(mu2) * np.sin(mu2)\n",
    "L_2 = (B - X) * np.cos(mu3) * np.sin(mu3)\n",
    "N = mu1 * np.cos(mu2) * np.cos(mu2)\n",
    "M = - (B - X) * np.cos(mu3) * np.cos(mu3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create mesh with gmsh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.initialize()\n",
    "gmsh.model.add(\"mesh\")\n",
    "p0 = gmsh.model.geo.addPoint(0.0, X, 0.0, lcar)\n",
    "p1 = gmsh.model.geo.addPoint(L - mu4, X, 0.0, lcar)\n",
    "p2 = gmsh.model.geo.addPoint(L, X, 0.0, lcar)\n",
    "p3 = gmsh.model.geo.addPoint(L + mu6 - L_2, H_2 + M, 0.0, lcar)\n",
    "p4 = gmsh.model.geo.addPoint(L + mu6, H_2, 0.0, lcar)\n",
    "p5 = gmsh.model.geo.addPoint(L, B, 0.0, lcar)\n",
    "p6 = gmsh.model.geo.addPoint(L + mu5, H_1, 0.0, lcar)\n",
    "p7 = gmsh.model.geo.addPoint(L + mu5 - L_1, H_1 + N, 0.0, lcar)\n",
    "p8 = gmsh.model.geo.addPoint(L, Y, 0.0, lcar)\n",
    "p9 = gmsh.model.geo.addPoint(L - mu4, Y, 0.0, lcar)\n",
    "p10 = gmsh.model.geo.addPoint(0.0, Y, 0.0, lcar)\n",
    "l0 = gmsh.model.geo.addLine(p0, p1)\n",
    "l1 = gmsh.model.geo.addLine(p1, p2)\n",
    "l2 = gmsh.model.geo.addLine(p2, p3)\n",
    "l3 = gmsh.model.geo.addLine(p3, p4)\n",
    "l4 = gmsh.model.geo.addLine(p4, p5)\n",
    "l5 = gmsh.model.geo.addLine(p5, p6)\n",
    "l6 = gmsh.model.geo.addLine(p6, p7)\n",
    "l7 = gmsh.model.geo.addLine(p7, p8)\n",
    "l8 = gmsh.model.geo.addLine(p8, p9)\n",
    "l9 = gmsh.model.geo.addLine(p9, p10)\n",
    "l10 = gmsh.model.geo.addLine(p10, p0)\n",
    "l11 = gmsh.model.geo.addLine(p1, p9)\n",
    "l12 = gmsh.model.geo.addLine(p2, p5)\n",
    "l13 = gmsh.model.geo.addLine(p5, p8)\n",
    "line_loop_rectangle_left = gmsh.model.geo.addCurveLoop([l0, l11, l9, l10])\n",
    "line_loop_rectangle_right = gmsh.model.geo.addCurveLoop([l1, l12, l13, l8, -l11])\n",
    "line_loop_bifurcation_top = gmsh.model.geo.addCurveLoop([l5, l6, l7, -l13])\n",
    "line_loop_bifurcation_bottom = gmsh.model.geo.addCurveLoop([l2, l3, l4, -l12])\n",
    "rectangle_left = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_left])\n",
    "rectangle_right = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_right])\n",
    "bifurcation_top = gmsh.model.geo.addPlaneSurface([line_loop_bifurcation_top])\n",
    "bifurcation_bottom = gmsh.model.geo.addPlaneSurface([line_loop_bifurcation_bottom])\n",
    "gmsh.model.geo.synchronize()\n",
    "gmsh.model.addPhysicalGroup(1, [l10], 1)\n",
    "gmsh.model.addPhysicalGroup(1, [l0, l1, l2, l4, l5, l7, l8, l9], 2)\n",
    "gmsh.model.addPhysicalGroup(1, [l3, l6], 3)\n",
    "gmsh.model.addPhysicalGroup(1, [l11], 4)\n",
    "gmsh.model.addPhysicalGroup(2, [rectangle_left], 1)\n",
    "gmsh.model.addPhysicalGroup(2, [rectangle_right], 2)\n",
    "gmsh.model.addPhysicalGroup(2, [bifurcation_top], 3)\n",
    "gmsh.model.addPhysicalGroup(2, [bifurcation_bottom], 4)\n",
    "gmsh.model.mesh.generate(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert to a dolfinx mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh, subdomains, boundaries = gmsh_to_fenicsx(gmsh.model, gdim=2)\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save final mesh, subdomains and boundaries (dolfinx format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with XDMFFile(MPI.COMM_WORLD, \"bifurcation.xdmf\", \"w\") as file:\n",
    "    file.write_mesh(mesh)\n",
    "    file.write_meshtags(subdomains)\n",
    "    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)\n",
    "    file.write_meshtags(boundaries)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
